Interaction-induced negative differential resistance in asymmetric molecular junctions 
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Combining insights from quantum chemistry calculations with master equations, we discuss a 
mechanism for negative differential resistance (NDR) in molecular junctions, operated in the regime 
of weak tunnel coupling. The NDR originates from an interplay of orbital spatial asymmetry 
and strong electron-electron interaction, which causes the molecule to become trapped in a non- 
conducting state above a voltage threshold. We show how the desired asymmetry can be selectively 
introduced in individual orbitals in e.g., OPE-type molecules by functionalization with a suitable 
side group, which is in linear conjugation to one end of the molecule and cross-conjugated to the 
other end. 



PACS numbers: 

I. INTRODUCTION 

The field of single-molecule electronics has been ex- 
panding rapidly during recent years, as techniques to 
electrically contact and control single molecules in a 
transport junction have improved [lH3- studying the 
electric current, /, through the molecule as function of 
the applied bias voltage, V, spectroscopic information 
can be extracted • In three-terminal setups a gate volt- 
age, Vg, can be used to control the electrostatic potential 
on the molecule, allowing for detailed spectroscopy . 

A long standing problem in molecular electronics is 
how to design molecules with reproducible and distinc- 
tive transport characteristics, e.g., diode-like rectifica- 
tion as in the original proposal of Aviram and Rat- 
ner @. Another interesting non-linear transport char- 
acteristic is negative differential resistance (NDR), i.e., a 
current which decreases with increasing voltage. NDR 
is known to occur in resonant tunneling diodes (RTDs) 
and has in this context been suggested for a large num- 
ber of device applications, such as high-frequency oscil- 
lators , logic circuits [1| , and analog-to-digital convert- 
ers 0. Also molecular junctions have been found to 
exhibit NDR mediated by a mechanism similar to the 
working principle of RTDs: a voltage-induced alignment 
of narrow bands in the electrode and molecule is followed 
by a mis-alignment at higher voltages, causing the cur- 
rent level to drop [loj . Such mechanisms can, however, 
not explain the huge NDR effect observed in different 
derivatives of oligo(phenyleneethynylene) (OPE), first 
observed by Chen et. al. Here explanations have 

instead involved either voltage-induced conformational 
changes [H, [11], or charge- injection [ll|, 113, Ell, both 
mechanisms forcing the molecule into a less conductive 



state above a certain voltage threshold. Clearly, under- 
standing the origin of NDR in e.g., such OPE molecules, 
being distinctly different from that of standard RTDs, is 
of great importance. The situation is, however, compli- 
cated by inconsistencies in the experimental findings, and 
while some experiments observe NDR in individual OPE- 
type molecules [l^ , others seem to indicate that NDR is 
only observed when large assemblies of such molecules 
are measured ■ 

In the above works the tunnel coupling between 
molecule and electrodes was rather strong. In gated 
(three-terminal) single- molecule junctions , the tun- 
nel coupling is often much smaller, on the same scale 
as the thermal energy ksT. NDR in this parameter 
regime has been theoretically discussed in various differ- 
ent types of molecules and can originate from e.g, Pauli 
spin-blockade [l^ (a well-known tool for spin-detection 
in semi-conducting quantum dot systems [Tol. [20| ) . vibra- 
tional blockade [2ll. [22| , or combinations of spin and vi- 
brational effects [23]. Additionally, interference between 
orbitally degenerate states has been predicted to give rise 
to NDR [2J] (analogous to quantum dots with ferromag- 
netic electrodes [2^ ) . There are also several experimental 
observations of NDR in molecular junctions operated in 
this regime, see e.g., Refs. [26l. [27|. 

In this work, we show that intrinsically asymmetric 
molecular junctions, such as OPE functionalized with ap- 
propriately chosen side groups, give rise to I{V) curves 
showing large NDR effects for low temperatures and weak 
tunnel coupling. Above a voltage threshold, a charge 
carrier (electron or hole) can tunnel onto the molecule 
and get trapped in a molecular orbital (MO) which has 
a significant overlap with the lead states of one con- 
tact only. The Coulomb interaction prevents additional 
charge transfer to the molecule, thereby blocking trans- 
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FIG. 1: Illustration of the NDR mechanism. Shown are 
both the molecular structure, with frontier LUMO (bottom) 
and LUMO+1 (top) orbitals, and a level-diagram, indicating 
the occupation of these orbitals. The HOMO (not shown) is 
filled with two electrons and transport involves the neutral 
and singly reduced molecular charge states. Rectangles in- 
dicate the source and drain conduction bands with chemical 
potentials fXs and fid- Below a (negative) threshold voltage, 
-|Vth|, an electron tunneling out to the drain from the LUMO 
can be followed by one tunneling into the LUMO -I- 1 from 
the source. Due to the orbital asymmetry, the electron in 
the LUMO-f 1 cannot leak out into the source and, due to 
the strong Coulomb interaction, it blocks additional charge 
transfer to the molecule, leading to a dramatically reduced 
current for V < — |Vth|. 



port and leading to NDR, sec Fig.[TJ Thus, our proposed 
mechanism arises because of strong Coulomb interaction 
between electrons on the molecule. Despite screening by 
conduction electrons in the leads, the Coulomb energy 
associated with charging molecules similar to those con- 
sidered here is indeed large (hundreds of meV per added 
electron or hole) [I,!!!]- Nonetheless, electron-electron 
interaction is often neglected, or treated on a mean- field 
level, in popular approaches to transport theory of molec- 
ular junctions, which would therefore fail to correctly de- 
scribe the effect discussed here. We instead use a master 
equation (ME) approach and calculate the nonequilib- 
rium occupations of the molecular many-body (i.e., in- 
cluding electron-electron interaction) energy-eigenstates, 
\Na), where is the number of valence electrons on 
the molecule and a labels the different iV-clcctron eigen- 
states. From these occupations, the current can be calcu- 
lated in the non-linear regime. Thus, both Coulomb in- 
teractions and the non-equilibrium condition are included 
non-perturbatively. With knowledge of the mechanism 
underlying the NDR, we use density functional theory 
(DPT) to investigate different OPE-type derivatives and 
find molecular structures where the effect is maximally 
pronounced. An ab initio approach to transport based 
fully on DPT is inadequate for our purposes because of 
the mean-field treatment of interactions. However, DPT 
can be used to predict the spatial structure of the (single 
particle) MOs, ^i. DPT can thus provide input param- 
eters to a many-body transport theory such as ME. In 
Appendix |^ we provide the details of the ME approach, 
as well as the link to the DPT calculations, i.e., the con- 
nection between -0^ and \Na). We also show how DPT 
can be used to estimate the the ratio jI/Jj, where •jI 



is the hybridization strength between MO -ipi and the 
continuum of electron states in electrode r, and we give 
the relation between 7^ and the many-body tunnel cou- 
plings, ^,^1 , determining the time-scale associated 
with a transition between states \Na) and \N'a') as a 
result of electron tunneling. 



II. ORBITALS AND STABILITY DIAGRAMS 

To identify molecules which are expected to show 
strong NDR effect, we have performed DPT calculations 
on a large number of OPE-type (and similar) molecules, 
three of which are discussed below. By choosing differ- 
ent side groups, asymmetry can be selectively introduced 
into different orbitals and thereby give rise to block- 
ing states and NDR at different gate and bias voltages. 
Electron-accepting substituent groups are expected to in- 
fluence most strongly the LUMO and LUMO-t-l, while 
electron-donating substituent groups are expected to in- 
fluence most strongly the HOMO and HOMO-1. MOs 
were calculated at the B3LYP/6-311G(2d,p) level using 
the Gaussian 03 Program package. Visualizations of the 
orbitals for a selected group of interesting molecules are 
shown in Pigs. [IHH 

Along with the orbital structure we also show the result 
of our ME calculations. Our purpose here is not to make 
precise predictions for a specific molecule, but rather to 
illustrate the general connection between NDR and the 
asymmetry of given orbitals. Rather than using the or- 
bital energies calculated with DPT (see Appendix [A)) , 
we therefore everywhere assumed energy-equidistant or- 
bitals, HOMO-1, HOMO, LUMO and LUMO-hl, with 
energy separation AE. All higher/lower orbitals are as- 
sumed empty /filled and are not included in the calcula- 
tions, and we focus on the singly oxidized, neutral, and 
singly reduced charge states only (for notational sim- 
plicity we always refer to the orbitals by their names 
in the neutral molecule). The number of valence elec- 
trons included in the ME calculations is thus = 3, 4, 5, 
and the many-body eigenstates are constructed by dis- 
tributing those electrons over the four orbitals as dis- 
cussed in Appendix |^ We assume the charging energies 
Uij between electrons in MO ipi and ipj all to be equal, 
Uij = 3AE. In all ME calculations we used the ther- 
mal energy fc^T = AE/50. Realistic values might be 
AE ~ 100 meV, which would mean T - 20 K. Typically, 
temperatures < 4 K can be reached in experimental se- 
tups, but we use a larger value to demonstrate the sta- 
bility of the NDR against thermally activated processes. 
Purthermore, in the regime of intermediate tunnel cou- 
pling, ksT < P, tunnel broadening (which is not included 
in the ME) would broaden the conductance features and 
reduce the NDR in a similar way as a larger temperature. 
Pinally, we simply take the electrode-hybridizations of 
the MOs to reflect the observed orbital asymmetry, i.e., 
for asymmetric orbitals they are taken much smaller at 
the lead attached to the side of the molecule with a small 
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orbital weight. The tunnel couplings of the many-body 
states are then calculated from Eq. (|A1|) in Appendix [K\ 

The results of the ME calculations are conveniently 
presented as so-called stability diagrams, i.e., as differ- 
ential conductance, dl/dV, plotted as function of gate 
and bias voltage, see e.g.. Fig. [DJb). Zero gate voltage 
has been chosen to correspond to the electrode Fermi 
levels (at zero bias) being halfway between the HOMO 
and LUMO of the neutral molecule (in reality this align- 
ment depends not only on the molecule, but also on 
the details of the electrodes and electric circuit). The 
left /right crossing point at zero bias corresponds to ox- 
idation/reduction of the molecule and the distance be- 
tween those points is the difference between the ioniza- 
tion energy and the electron affinity. The linear (small 
bias) conductance thus give similar information as can be 
obtained from voltammetry. The scaling factor a (gate 
coupling) depends on the capacitances of the tunnel junc- 
tions [201 . The gate-dependent lines at finite bias corre- 
spond to the energy condition that an (additional) tunnel 
process becomes possible. This happens when a molecu- 
lar chemical potential, A*7Va,(7V-i)Q' = Epfa — £'(Ar-i)Q', 
where En a is the energy of the state \Na), comes into the 
bias window. Usually, this leads to a stepwise increase of 
the current, giving a peak in dl/dV. From the voltage 
positions of these lines one can thus extract the many- 
body excitation energies and a three-terminal molecular 
junction therefore works as a spectroscopic tool. The 
slopes of the lines depend on the capacitances of the tun- 
nel junctions, which can therefore be extracted in an ex- 
perimental situation and used to calculate the gate cou- 
pling a poj . For simplicity we here assume the capaci- 
tances associated with the source and drain contacts to 
be equal, and equal for all states, and furthermore ap- 
ply the bias voltage in a symmetric fashion, fis = 
fid = —eV/2^ where — e is the electron charge, leading to 
a uniform slope of all dl/dV lines. 

In the first example, shown in Fig.[5][a), OPE has been 
functionalized by a cyano (CN) group. The weights of 
all orbitals are almost left-right symmetric, with the ex- 
ception of the LUMO-l-1, which almost completely lacks 
weight on one of the end benzene rings, namely the end 
which is in linear conjugation to the cyano group (ortho 
configuration). The other end, with all the LUMO-l-1 
weight, is cross-conjugated to the cyano group (meta 
configuration). This dependence of the LUMO and 
LUMO-l-1 on the conjugation pathway to the substituent 
group is also evident from other structures described be- 
low. In general, the LUMO always has weight at the " or- 
tho-end^\ while the LUMO-l-1 always has weight on the 
" meta-cnd" of the molecule. Figure [2][b) shows the sta- 
bility diagram obtained from the model described above. 
For positive gate and negative bias voltage, there is a 
clear negative (blue) dl/dV line, i.e, NDR. The red solid 
line in Fig. [2jc) shows the current as function of bias 
voltage along the green dashed line in (b). At this value 
of Vg, for very small \V\^ the molecule is reduced and 
the current suppressed since interaction effects (Coulomb 
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FIG. 2: (a) Molecular structure and frontier orbitals for OPE 
with a CN side group, (b) Stability diagram, i.e., dl /dV as 
function of Vg and V . The DFT calculations predict a very 
asymmetric LUMO-l-1 and we set all electrode hybridizations 
equal to 7, except 7lumo+i = 7/10- (c) Red solid curve 
shows source-drain electron current, I{V), along the dashed 
green line in (b). Green dashed and blue dotted curves show 
the same, but for a larger coupling of the LUMO-f 1 to the 
drain, 7/4 and 7/2 respectively. 



blockade) and the discrete nature of the molecular or- 
bitals prevent charging/discharging of the molecule. In 
the regime of weak tunnel coupling discussed here, the 
molecular orbitals are very narrow in energy and trans- 
port including only virtual charging of the molecule is 
suppressed. When \V\ becomes larger, the molecule can 
fiuctuate between the neutral and reduced charge states 
as electrons tunnel into and out of the LUMO, and a fi- 
nite current starts to flow. For even more negative U, the 
current then goes down (NDR) as sketched and described 
in Fig. [T] when the molecule can be reduced by an elec- 
tron tunneling from the drain into the LUMO-l-1, rather 
than LUMO, this electron is "trapped" since it cannot 
tunnel out into the source due to the orbital asymmetry; 
additionally, it blocks transport through the other or- 
bitals, since the strong Coulomb repulsion prevents more 
electrons from tunneling onto the molecule. The block- 
ade is lifted and the current recovers when the voltage 
allows an electron to tunnel out to the source from the 
HOMO instead. For F > there is no NDR effect and 
the LUMO+l is simply not visible in the transport spec- 
trum. The current in the blockaded region recovers if the 
LUMO-l-1 hybridizes more strongly with the drain elec- 
tron states, as is evidenced by the green dashed and blue 
dotted curves in Fig.[21[c). Note that the three curves are 
almost identical at positive bias. 

Next, in Fig. [H^a), we show OPE functionalized by a 
dicyanoethylene group, which leads to a strongly asym- 
metric LUMO. Analogous to the above case of an asym- 
metric LUMO-l-1, this gives rise to a pronounced NDR, 
see Fig. [31^b). In this case, however, the blocking state is 
an excited state of the neutral molecule, where an clcc- 
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FIG. 3: (a) Molecular structure and frontier orbitals for OPE 
with a dicyanoethylene side group, (b) Stability diagram, 
calculated using all lead hybridizations equal to 7, except 
7lumo — 7/10- (c) Red solid and green dashed curves show 
I{V) along the corresponding lines in (b). 



tron has been promoted from the HOMO to the LUMO, 
and therefore the NDR is seen instead at negative gate 
voltages (that it occurs at positive bias is simply because 
there is a small coupling of the LUMO to the drain, rather 
than source; in an experimental situation this is not easily 
controllable since it depends on how the molecule binds 
inside the junction). The red solid curve in Fig. [3jc) 
shows the current along the red solid line in (b). There is 
no strong NDR at positive gate voltage. Instead, the low- 
bias current is strongly suppressed, the reason being that 
the asymmetric LUMO is occupied in the ground state of 
the reduced molecule, which therefore acts as a blocking 
state. Here the I{V) characteristic instead look similar 
to a regular diode, being strongly asymmetric between 
positive and negative bias, see the green dashed curve in 
Fig-EJc) which shows the current along the green dashed 
line in (b). Interestingly, using the gate voltage, one can 
control whether the positive or negative bias direction 
gives rise to a large current. 

Asymmetric HOMO or HOMO-l orbitals would give 
similar stability diagrams as above, but with NDR due 
to blocking hole, rather than electron, states. 

In OPE-type molecules, the central ring can rotate 
around the triple bonds and the ~ 90° rotated species 
have only a slightly higher ground state energy compared 
to the planar ones. Ref. [Ijl, which studied molecules 
similar (although not identical) to those discussed above, 
concluded that such rotations can change the spatial 
structure of the orbitals and a conformational change 
can switch the molecule between conducting and non- 
conducting states. A similar explanation was suggested 
in the experimental work of Ref. [l^. However, it is not 
at all clear that such rotations can take place once the 
molecule resides inside the junction, possibly "lying" on 
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FIG. 4: (a) Molecular structure and frontier orbitals for a 
rigid molecule, where the benzene rings are fused to five- 
membered rings, functionalized with a NO2 group, (b) Sta- 
bility diagram, calculated using all lead hybridizations equal 
to 7, except 7LUMO = 7homo-i = 7/10. 



top of the back gate. We emphasize that all our consider- 
ations are for the lowest energy geometry, which is (close 
to) planar, and the predicted NDR effect does not rely on 
conformational changes. Indeed, we are able to create a 
completely rigid molecule, in which three benzene rings 
are "locked" in a planar conformation by fusing them 
to two five-membered rings (corresponding to two fluo- 
rene sub-structures). This molecule shows spatial orbital 
asymmetry when the central ring is functionalized with 
a NO2 group (Fig. 4a). As in Fig. [3Ib), the asymmetric 
LUMO suppresses transport close to the right crossing 
point and induces NDR at positive bias close to the left 
one, see Fig.H^b). Additionally, the asymmetric HOMO- 
1 gives rise to NDR at negative bias, as tunneling by a 
hole into this orbital blocks hole-transport through the 
molecule. 



III. CONDITIONS FOR NDR 

The type of NDR effect discussed above relics on a 
meta-stable (blocking) state being occupied. Therefore, 
if this state can relax by other means than tunneling, e.g., 
by photon emission, the degree of blocking and therefore 
the magnitude of the NDR will be diminished. The life- 
time should here be compared with the typical time be- 
tween electron tunneling events in the unblocked regime, 
i.e., with e/I. Even for weakly coupled molecules, this is 
a rather short time, e.g., a current / nA, corresponds 
to a time-scale ^ 10^^° s, which should indeed exceed 
the lifetime of most excited states. That excited molecu- 
lar states are indeed stable on such time-scales was also 
the conclusion of the experimental work in Ref. [l^] . 

Our conclusions concerning the orbital asymmetry 
were based on the chemistry of the isolated molecule, 
i.e., neglecting the influence of the bonds to the elec- 
trodes, the details of which depend on the (unknown and 
uncontrollable) exact bonding geometry. However, we 
here consider the weak coupling regime, meaning that the 
molecule-electrode interactions arc due to physisorption. 
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rather than chemisorption, and the resulting perturba- 
tions of the molecular orbital structure should therefore 
be small. Furthermore, as is seen in Fig. [2ljc) , it is suffi- 
cient for NDR if one coupling amplitude is twice as large 
as the other, so the effect should persist even if the in- 
trinsic asymmetry is partially destroyed by the binding 
to the electrodes. 

In addition, the proposed mechanism relies on the ends 
of the molecular chain binding to the electrodes, rather 
than some other configuration. This can at present not 
be experimentally controlled and therefore perhaps only 
a subset of all junctions incorporating such molecules ac- 
tually show NDR. It may be advantageous to function- 
alize the molecules with appropriate end-groups. If the 
electrodes are made of gold, one might e.g., use thiol end- 
groups, connected to the molecular chain through satu- 
rated carbon atoms, as in a methylene spacer, to keep 
the effective molecule-electrode coupling small (29j . 



IV. SUGGESTION FOR AB INITIO 
CALCULATIONS 



Finally, we suggest the possibility to investigate trans- 
port through the molecules studied here using a com- 
pletely ab initio approach. Due to the intrinsic many- 
body nature of the predicted NDR effect, the standard 
ab initio description of transport, based on combining a 
mean-field description of the molecule (e.g., using DFT) 
with non-equilibrium Greens functions (NEGFs), can- 
not be expected to yield correct results. Nonetheless, 
it might be of interest to undertake such calculations 
for comparison as they may well capture the existence 
of NDR. Within a mean-field theory, the Coulomb in- 
teraction will increase the energy, Ej , of the conducting 
orbital, , once the blocking orbital, ipi , enters into the 



bias window: 



Ej (X Uij{ni), 



where Ua is the Coulomb 



interaction between electrons in orbitals "04 and tpj [cf., 
Eq. (|A8|) in Appendix E] . and {rii) is the average occu- 
pation of orbital ipi. This will reduce the current carried 
by orbital ipj and might lead to NDR. Important here is 
that the non- equilibrium aspects of transport are treated 
correctly, i.e., that the transmission function, T{E,V), 
is calculated for finite V based on the non-equilibrium 
charge distribution, as is done e.g., in the implementa- 
tions discussed in Refs. [33,[3l|. Furthermore, the correct 
description of transport through the types of molecules 
discussed here, and their intrinsic NDR effect, would pro- 
vide an interesting testing ground for newly developed ab 
initio approaches which attempt to go beyond a mean- 
field description, see e.g., Refs. j33 - l34| . 



V. SUMMARY AND CONCLUSIONS 

We have shown that an interplay of orbital spatial 
asymmetry and strong electron-electron interaction can 
lead to NDR in molecular junctions. Using DFT we 



showed that functionalization of OPE-type molecules 
with suitable side groups introduces the desired asym- 
metry in selected individual orbitals. Using the chemical 
insights as a starting point, a master equation approach 
was used to predict the resulting three-terminal transport 
signal. 

Although the purpose of this paper was to investigate 
the transport behavior in the regime of weak tunnel cou- 
pling, we mention the possibility that our results may still 
explain some of the NDR effects seen in similar molecules, 
e.g., in Refs. [O, [13, EBl , despite the larger tunnel cou- 
pling in those experiments. However, we emphasize that 
NDR has also been seen in OPE-type molecules with- 
out orbital asymmetry. Experimental investigation of the 
type of molecules consider here in three-terminal junc- 
tions would be very interesting, as the increased level of 
control achieved with the gate electrode would allow more 
careful comparison to our predictions, and additionally 
provide the possibility of tuning the device to the regime 
of strong NDR. 
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Appendix A: Master equations for the molecular 
junction 

We consider a molecule which is coupled to a source 
(s) and a drain (d) electrode. Neglecting at first the 
molecule-electrode hybridization, the molecule is char- 
acterized by a set of many-body energy eigenstates \Na), 
where N is the number of valence electrons on the 
molecule (the number of electrons wc explicitly include 
in the calculations) and a labels the different A^-clcctron 
eigenstates (a here labels also the spin degrees of free- 
dom). We let {tpi} be a single particle orbital basis for 
the molecule and define operators cl^ which create an 
electron with spin-projection a in orbital ipi (if needed, 
one could also include an index a on the orbitals, allow- 
ing e.g., for different orbitals for a and /3 spin in spin 
unrestricted calculations). Furthermore, electrons in or- 
bital ipi hybridize with the continuum of electron states 
in electrode r = s,d with amplitude ^l^. The tunnel rate 
between molecular eigenstates \Na) and \N'a') is then 
defined as 



Na,N'c 



— dr 

n 



\N'a') 



(Al) 
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where only the first (second) term gives a non-zero con- 
tribution it N' ^ N -1 {N' = N + 1). Here 4 is the 
density of states in electrode r, which we assume to be 
energy-independent (an energy dependence is not needed 
for the NDR mechanism discussed here and would change 
the results only quantitatively). 

We consider the regime of weak molecule-electrode 
tunnel coupling. Here "weak coupling" means that all 
hr are small compared to the energy separation of the 
states involved and, in addition, smaller than the thermal 
energy scale fc^T. In this regime, stationary state (time- 
independent) transport can be described by leading or- 
der perturbation theory in F, yielding rate (or master) 
equations [11] , which are solved to obtain the occupation 
probabilities P/vq of the molecular states: 

a.,N'a 'Pn'c' — Wn'c' ,NaPNa) , (A2) 

N'a' 

1 - ^Ptvc, (A3) 

Na 

where Eq. (|A3p enforces probability normalization. 
Equation (|A2[) simply states that the occupation of a 
state is determined by the balance of in-going and out- 
going processes (the sum of all tunnel processes going 
into state |A^a), weighted by the occupation of the corre- 
sponding initial state \N'a'), minus all tunnel processes 
going out of state \Na)). In leading order the rate ma- 
trix W oc r describes only tunnel processes involving a 
single electron. It thus has non-zero elements only when 
N and N' differ by exactly one, and those elements are 
given by 

WNa,(N-l)a' = '^^^'n a.{N -\)a' 

r 

X friENa - -E(w-l)a'), (A4) 

WNaXN+l)a' = ^^NaAN+Da' 
r 

X [l - fr{E^N+l)a' - ENa)] , (A5) 

where fr{E) = l/(e(^"^'-)/'=B^ -f- 1) is the Fermi function 
of electrode r with chemical potential /i^ and Ej^a is the 
energy of the molecular eigenstate \Na). 

Solving Eqs. (|X2|) - (|X3)) with the rates (|X4|) - ((M)) gives 
the non-equilibrium occupation probabilities PNa of all 
molecular many-body eigenstates. Based on these occu- 
pations, the current through the molecular junction can 
be calculated. The source-drain electron current can be 
obtained from 

^ = Y.J2^Na.N'a'PN'a', (A6) 

Na N'a' 

where the current rate matrix W-^ is similar to (|A4p - 
(|A5p . but including only tunnel processes involving elec- 
trode r ~ s and a plus/minus sign for electrons tunneling 
onto/out of the molecule. 

The needed input to the transport calculations is thus 
the full spectrum of molecular many-body eigenstates. 



\Na), with energies Enc, as well as the hybridizations 
7fg. between the single-particle orbitals ipi and the lead 
electron states. The many-body states can be expressed 
as a sum of single-particle states (Slater determinants in 
first quantization language) with N electrons: 

\Na)= Cl;r-^ct,---4.j0), (A7) 

where the sum runs over all collections {ik<^k}^=i satis- 
fying ii < ■ ■ ■ < iN and |0) is the state with zero valence 
electrons. From the hybridizations 7^^, we can thus cal- 
culate the tunnel couplings by inserting (|A7|) into (jAip . 
However, finding the eigenstates (equivalent to finding all 

^Na'^^^)^ in principle involves diagonalizing the many- 
body Hamiltonian of the molecule, excluding tunnel cou- 
pling to the electrodes, but including the electrostatic 
environment of the electrodes, typically leading to re- 
duced Coulomb repulsion and HOMO-LUMO gap due to 
screening by conduction electrons in the electrodes [28t . 
This is clearly a very difficult problem. In the follow- 
ing wc instead suggest a simple approximative approach, 
which still combines insights from ab initio quantum 
chemistry calculations with non-pcrturbative treatment 
of electron-electron interactions, crucial to correctly de- 
scribe transport in the weak tunneling regime. 

The basis states ipi are taken as the Kohn-Sham (KS) 
orbitals of the neutral molecule, calculated using DFT, 
and thus already include effects of electron-electron in- 
teractions on a mean-field level. We now assume that we 
can construct each many-body state by simply distribut- 
ing the N valence electrons over these orbitals. i.e., each 
state is given by setting one c\^^''^ = 1 and all others 
to zero in (jATp . This is justified when the MOs are not 
significantly changed upon charging, which is often the 
case. However, as long as the orbital asymmetries are 
not completely destroyed upon charging, the predicted 
NDR effect is not quantitatively affected. (If, on the 
other hand, the change of the MOs as the molecule is 
charged is crucial, the above approach can simply be re- 
placed by one where each state \Na) is constructed in- 
stead from orbitals ijjNi calculated for the appropriately 
charged molecule.) The "many-body" states are thus ef- 
fectively single-particle states given by a single Slater de- 
terminant. The non-trivial effect of the electron-electron 
interaction instead enters in the corresponding energies 
ENa 1 which arc approximated by the sum of the energies 
Cicr of the occupied orbitals (the collection of which we 
call ONa), plus a Coulomb charging term: 

ENa = "i'^^^ H (^^) 

where Uij is the Coulomb charging energy associated 
with one electron in orbital i and one in orbital j. Uij 
could be calculated within various approximations. How- 
ever, the exact value depends on microscopic details of 
the junction [28| (not only the molecule), which are at 
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present uncontrollable, and unknown, in experimental 
situations. Fortunately, even though including a large 
Coulomb charging energy is crucial for a correct descrip- 
tion of transport, the exact value typically has only a 
quantitative effect and can be seen as a fitting parameter 
when comparing with experimental data. 

With the above approximations, the tunnel cou- 
plings dSH) simplify to r'-„_^,„, = 27rdr\-flJ^ if \Na) is 
obtained from | A'^'a') by adding or removing one electron 
with spin-projection a in orbital i, and = if 

\Na) and \N'a') do not differ only by the occupation of a 
single orbital. The remaining problem is then to find the 
hybridization parameters "fl^. Since tunnel amplitudes 
decrease exponentially with increasing barrier, it is rea- 
sonable to assume that tunneling only takes place into 
one or a few atoms at each end of the molecule which 
form the bond to the electrode. The natural basis is 
therefore the atomic orbitals (AOs), (j>j. We express the 
MOs in terms of the AOs, tpi = J2j '^)'t>ji which finally 
gives the tunnel coupling 

T^^^.N'o.'^'^^dr\Y^a)%,\\ (A9) 

j 



were 7jj^ is the hybridization of AO j with electrode r and 
\Na) and \N' a') differ by an electron in orbital i. In prin- 
ciple, one could attempt to calculate using ab initio 
methods. However, the result would again strongly de- 
pend on the (unknown and uncontrollable) details of the 
junction. In addition, in the weak coupling regime con- 
sidered here, only the ratio of the source and drain tunnel 
couplings affect the transport behavior, while their over- 
all magnitude only set the scale for current. The coeffi- 
cients = {ipi\4'j) can simply be read out of a DFT (or 
other quantum chemistry) program, provided that local- 
ized atomic orbitals are used as a basis. We note that if 
the KS orbitals would change substantially upon charg- 
ing, such that we would be forced to use iV-dependent 
orbitals, ^Ni, orbitals in different redox states are not 
orthogonal and Eq. (jA9p should be modified to contain 
a sum over all orbitals, weighted by the corresponding 
overlaps. 
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